% Generate a number of networks for the Machine Learning training dataset

%% Settings
% % Dataset
% % Number of training networks to generate
Ntr=2;
Nimg=1;
% rMeanBounds=[20 70]
% rStdBounds=[0 10]
% thMeanBounds=[5 20]
% thStdsBounds=[0 5]

% Network
poreNetInfo.DomainW=1000;
poreNetInfo.DomainH=1000;

poreNetInfo.nPor=100;
poreNetInfo.porDist='normal';
poreNetInfo.thDist='normal';

meanRPor=30;
devRPor=5;

poreNetInfo.porShape='circle';
poreNetInfo.thShape='smooth';

poreNetInfo.thDistParam.meanTh=10;
poreNetInfo.thDistParam.devTh=1;

poreNetInfo.MergeOverlappingPores=true;

frame=ceil(meanRPor+3*devRPor);

poreNetInfo.porDistParam.meanDPor=meanRPor*2;
poreNetInfo.porDistParam.devDPor=devRPor*2;

%% Generation of training samples
% Preallocate
poreData=cell(Ntr,1);
thData=cell(Ntr,1);


for nn=1:Ntr
    warning off;
    
    % Random seed control for repeatability
    rng(1);
    
    % Create parallel pool if none exists.
    gcp;
    
    tic
    
    bAlpha=poreNetInfo.DomainW/40; % alpha radius of the network boundary
    
    % Generate Voronoi Tesselation for the primary porosity
    fprintf('Generating Voronoi tesselation..');
    
    % Generation of the random point cloud on which the network will be based
    poreNetInfo.nPor=round( (poreNetInfo.nPor + 33)/1.9 );
    prePor(:,1)=frame+(poreNetInfo.DomainW-2*frame)*rand(poreNetInfo.nPor,1);
    prePor(:,2)=frame+(poreNetInfo.DomainH-2*frame)*rand(poreNetInfo.nPor,1);
    
    % Generate the Delaunay and Voronoi tesselations
    TR=delaunayTriangulation(prePor);
    [centers,Regions]=voronoiDiagram(TR);
    
    % Remove out-of-domain regions
    centers(centers(:,1)>poreNetInfo.DomainW,:)=inf;
    centers(centers(:,2)>poreNetInfo.DomainH,:)=inf;
    centers(any([centers(:,1)<0, centers(:,2)<0],2),:)=inf;
    
    ii=1;
    needed=true;
    while needed
        if ismember(inf,centers(Regions{ii},:))
            Regions(ii)=[];
            ii=ii-1;
        end
        
        ii=ii+1;
        
        if ii>size(Regions,1)
            needed=false;
        end
    end
    
    % Exlude out-of-domain pores and update the 'Regions'.
    % Mask of in-domain pores
    maskRef=ismember(1:length(centers),cell2vec(Regions));
    
    % Refine region information
    cum=zeros(length(centers),1);
    cum(1)=maskRef(1);
    for ii=2:length(maskRef)
        cum(ii)=cum(ii-1)+maskRef(ii);
    end
    
    [newcode,oldcode]=unique(cum);
    
    newcode=newcode(2:end);
    oldcode=oldcode(2:end);
    
    nReg=length(Regions);
    
    for ii=1:nReg
        Regions{ii}=changem(Regions{ii},newcode,oldcode);
    end
    
    % Refine pore locations
    centers=centers(maskRef,:);
    
    % Construct edges from region information
    % Initialize the edges matrix
    % Count the anticipated number of edges
    numEdges=0;
    for ii=1:nReg
        numEdges=numEdges+numel(Regions{ii});
    end
    
    edges=zeros(numEdges,2);
    regID=zeros(numEdges,1);
    
    oo=1;
    
    for ii=1:nReg
        
        numVertInReg=numel(Regions{ii});
        
        for kk=1:numVertInReg
            if kk==numVertInReg
                edges(oo,:)=[Regions{ii}(kk) Regions{ii}(1)];
                regID(oo)=ii;
                oo=oo+1;
            else
                edges(oo,:)=[Regions{ii}(kk) Regions{ii}(kk+1)];
                regID(oo)=ii;
                oo=oo+1;
            end
        end
    end
    
    
    % Remove duplicate edges
    for ii=1:numEdges
        if ~isnan(edges(ii,1))
            mask1=all(edges==edges(ii,:),2);
            mask2=all(edges==fliplr(edges(ii,:)),2);
            mask12=any([mask1,mask2],2);
            mask12(ii)=false;
            edges(mask12,:)=NaN;
            edges(mask12,:)=NaN;
        else
            continue
        end
    end
    
    edges=edges(~isnan(edges(:,1)),:);
    
    % Coordinates of starting and ending points of each throat
    startsMidP=centers(edges(:,1),:);
    endsMidP=centers(edges(:,2),:);
    
    % Calculate total area of the domain (convex hull)
    outBound=convhull(centers(:,1),centers(:,2));
    
    totalArea=area(polyshape(centers(outBound,1),centers(outBound,2)));
    
    poreNetInfo.nPor=size(centers(:,1),1);   % Find number of pores
    poreNetInfo.nTh=size(startsMidP,1);   % Find number of throats
    
    fprintf('Done\n');
    
    %% Generate the sizes of pores and throats
    radii=abs(normrnd(meanRPor,devRPor,poreNetInfo.nPor,1));
    widths=abs(normrnd(poreNetInfo.thDistParam.meanTh,poreNetInfo.thDistParam.devTh,poreNetInfo.nTh,1));
    
    %% Merge the overlapping primary pores if requested by user
    if poreNetInfo.MergeOverlappingPores
        fprintf('Merging Overlapping Pores..')
        [centers,radii,startsMidP,endsMidP,widths,Regions]=PoreMerger2(centers,radii,startsMidP,endsMidP,widths,Regions,edges);
        fprintf('Done \n')
        % Recalculate the number of pores and throats
        poreNetInfo.nPor=size(centers,1);
        poreNetInfo.nTh=size(startsMidP,1);
    end
    
    % Delete the pores removed during merging (denoted by NaN)
    NaNMask=isnan(radii);
    radii=radii(~NaNMask);
    centers=centers(~NaNMask,:);
    poreNetInfo.nPor=length(radii);
    poreNetInfo.nTh=length(widths);
    
    % Inlet/Outlet flagging
    maskInlet=centers(:,1)<0.1*poreNetInfo.DomainW;
    maskOutlet=centers(:,1)>0.9*poreNetInfo.DomainW;
    
    
    %% Record the data into a delimited text file
    poreData=[centers, radii, maskInlet,maskOutlet];
    fileID = fopen('trainingPoreData.txt','w');
    fprintf(fileID,'X','Y','R','InletFlag','outletFlag');
    fprintf(fileID,'%6.2f %12.8f\n',netData);
    fprintf(fileID,'endOfInstance');
    fclose(fileID);

    
    %% Image generation
    if nn<=Nimg
        poreNetInfo.fileName=strcat('train',num2str(nn),'.dxf');
        
        % Generate pores
        fprintf('Generating pores..')
        % Primary
        switch poreNetInfo.porShape
            
            case 'circle'
                % Generate the circles
                x = circleGen(radii);
                
                % Plot the pores
                for ii=1:poreNetInfo.nPor
                    % Translate  the center of the circle to the given pore location
                    x{ii} = x{ii} + centers(ii,:);
                end
                
            case 'triangle'
                % Generate the squares
                x = triGen(radii);
                
                % Plot the pores
                for ii=1:poreNetInfo.nPor
                    % Translate  the center of the circle to the given pore location
                    x{ii} = x{ii} + centers(ii,:);
                end
                
            case 'square'
                % Generate the squares
                x = sqGen(radii);
                
                % Plot the pores
                for ii=1:poreNetInfo.nPor
                    % Translate  the center of the circle to the given pore location
                    x{ii} = x{ii} + centers(ii,:);
                end
                
        end
        fprintf('Done \n')
        
        % Generate throats
        fprintf('Generating throats..')
        
        % Calculate throat direction vectors
        throatVecs=[endsMidP(:,1)-startsMidP(:,1), endsMidP(:,2)-startsMidP(:,2)];
        throatLengths=(throatVecs(:,1).^2+throatVecs(:,2).^2).^0.5;
        throatVecs=throatVecs./repmat(throatLengths,1,2);
        
        switch poreNetInfo.thShape
            
            case 'smooth'
                % Generate rough rectangles
                xt = recGen(throatLengths, widths,throatVecs);
                
                % Translate  left-middle point of the rectangle to the desired location
                for ii=1:poreNetInfo.nTh
                    xt{ii} = xt{ii}+repmat( throatVecs(ii,:)*throatLengths(ii)/2+startsMidP(ii,:),size(xt{ii},1),1 );
                end
                
        end
        
        fprintf('Done \n')
        
        x=[x,xt];
        
        % Union the pores and throats
        fprintf('Performing polygon boolean operations..')
        [shapeCell , porArea] = unionMultiShapes(x, 512, 256);
        fprintf('Done \n')
        
        poreNetInfo.porosity = porArea / totalArea;
        
        %% Generate DXF file
        fprintf('Generating the DXF file..')
        genDXF(shapeCell, poreNetInfo.DomainW, poreNetInfo, poreNetInfo.fileName);
        fprintf('Done \n \n')
        fprintf('Total generation time: %1.2f seconds \n', toc)
    end
end